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ABSTRACT 

A finite difference formulation is presented for wave propagation in 
a rectangular two-dimensional duct without steady flow. The difference 
technique, which should be useful in the study of acoustically treated inlet 
and exhausts ducts used in turbofan engines, can readily handle acoustical 
flow field complications such as axial variations in wall impedance and 
cross section area. In the numerical analysis, the continuous acoustic 
field is lumped into a series of grid points in which the pressure and veloc 
ity at each grid point are separated into real and imaginary terms. An 
example calculation is also presented for the sound attenuation in a two- 
dimensional straight soft- walled suppressor. 

INTRODUCTION 

Because of increasingly stricter environmental noise regulations' 
noise suppressors are now an integral part of turbojet engine design and 
other industrial turbomachinery applications. In the particular example 
of the aircraft industry, turbojet suppressors must be optimized for mini- 
mum weight and length to decrease the direct operational cost of the air- 
craft and to increase the return on investment. 



Generally, aircraft suppressors dissipate acoustical energy by viscous 
action in walls containing absorbing materials. Presently, estimates for 
the attenuation of acoustic energy (noise) in suppressor ducts are based 
on the analytical study (refs. 1-7) of acoustic waves propagating into 
infinitely long straight ducts with uniform impedance (absorbers) along the 
walls. The analytical techniques, however, are limited to relatively 
simple geometries. 

At the present time, there is a need for more flexible suppressor 
design techniques which can handle acoustical flow field complications in 
engine ducts such as axial variations in wall impedance and cross sectional 
area as would occur in a sonic inlet. To meet the need, the present paper 
develops a numerical finite difference technique which can be used .in the 
prediction of sound attenuation in turbojet engine ducts as well as other 
environmental noise abatement problems. 

In particular the present paper formulates a solution to the Helmholtz 
equation with finite wall impedance by a difference technique. The differ- 
ence solutions bypass the conventional eigenvalue problem with its assoc- 
iated modes which have been considered in earlier works on noise propaga- 
tion in ducts (refs. 1-7). As a result of the difference formulation, the 
propagation of noise is treated as a diffusion process analogous to problems 
in thermodynamics involving heat flow or fluid dynamics involving the 
transport of vorticity. 

Immediately following the mathematical development of the difference 
technique, example solutions are presented for a comparison between the 
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finite difference and the analytical results for sound propagation in a one- 
dimensional hard wall duct and a two-dimensional soft walled duct without 
steady flow. 
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SYMBOLS 

speed of sound 
decibels, Eq. (19) 
acoustic power, Eq. (18) 
frequency 

channel height, see Fig. 1 
acoustic intensity 


L length of duct 

m total number of grid rows 

n total number of grid column 

P* pressure fluctuation, P’(x T , y', t) 

p dimensionless pressure fluctuation, p(x, y) = P'/P^ 

p^ amplitude of pressure fluctuation at duct entrance 

t time 

u dimensionless acoustical particle velocity 

x, y dimension axial coordinate, x'/H, and transverse coordinate, y T /H 

Ax, Ay grid spacing 
Z acoustic impedance 

^ specific acoustic impedance 

dimensionless frequency (eq. (4)) 

0 specific acoustic resistance 

k acoustic conductance ratio 
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X wavelength 

p density 

0 acoustic susceptance ratio 

X specific acoustic reactance 

co circular frequency 

Subscripts: 

’ prime indicates a dimensional quantity 

(1) real part 

(2) imaginary part 
Subscripts: 

e exit condition 

i, j i axial index, j transverse index, see Fig. 2 

t transverse y direction 

x axial position 

y transverse position 

GOVERNING EQUATIONS AND BOUNDARY CONDITIONS 

The operation of a typical turbojet suppressor using a Helmholtz 
resonator as a noise absorber is shown schematically in Fig. 1. The 
noise source shown in Fig. 1 has associated with it some initial pressure 
profile and some entering acoustic energy flux represented by the symbol 

1 (intensity). The entering acoustic energy is dissipated along the length 
of the duct by viscous action in the Helmholtz absorbers imbedded in the 
walls of the duct. The acoustic particle fluid velocities u’ and uV. shown 

x y 

in Fig. 1 play an important role in the dissipation process as well as in the 
transmission of the acoustic energy along the duct. The governing equations 
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and initial and boundary conditions which describe the transmission and 
absorption process will now be presented. 

WAVE EQUATION 

The governing differential equation describing the propagation of 
acoustical energy in the wave equation. The two dimensional form of 
the wave equation without steady flow is given by 

a 2 pr a 2 pf ! a 2 p , 

+ = ( 1 ) 

a ,2 . .2 2 .,2 

3x f dy ’ c dt 

In the above equation, the prime, ', is used to denote a dimensional 
quantity. These and all other symbols used in this report are defined 
in the list of symbols. As customary for steady state, the solution is 
assumed to be of the form 

P f (x\ y T , t) = p' (x r , y’) e +iwt ( 2 ) 

Substituting Eq. (2) into Eq. (1), and introducing the dimensionless 
parameters, p, x, y and rj yields the classic Helmholtz equation. 

^ ^ + (27tt?) 2 p = 0 (3) 

0 x 2 3y 2 

where the dimensionless frequency 77 is given as 

= J_ ^ = Hf = H 
2 t c c \ 


( 4 ) 
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The frequency parameter ij represents the ratio of duct height H to 
acoustic wavelength for the geometry shown in Fig. 1. The dimension- 
less height y ranges from zero to 1 while the dimensionless length x 
ranges between zero and L/H. 

When using the exponential notation displayed in Eq. (2), the dimen- 
sionless pressure p has in general both real and imaginary parts. Thus, 

P (x, y) = P^ (x, y) + ip^ (x, y) (5) 

Consequently, Eq. (3) can now be broken into its real and imaginary 
parts by substituting Eq. (5) into Eq. (3), 

(1) g2 (1) 2 (1) 

^-2 — + ^-2 — + (2tt77) V 1} = 0 (6) 

- 2 a 2 

dx By 


aV 2 > , aV 2 > 

ax 2 ay 2 


+ (27n7) 2 p (2) 



(V 


Equations (6) and (7) represent the basic governing equations for 
noise propagation which will be solved later. 

Acoustic Particle Velocity 

The boundary conditions in acoustics are generally given in terms 
of impedances which relate the pressure and velocity fields at the boundary 
(see schematic insert in fig. 1). For a harmonic solution of the form 
given by Eq. (2) and in the absence of convective velocities, the dimension- 
less momentum equations yield the following expressions for the scalar 
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velocities in terms of pressure gradients: 


. 3p . 3p 

u = 1 . — “ u = 1 — — 


3x 


3y 


where the dimensionless velocities are defined as 


pwH . pw H , 

u x- n u x u y " u y 

PA pa 


(8) 


(9) 


Impedance Boundary Condition 
The specific acoustic impedance ^ is defined as 

V p' 

= = at y = 0 and 1 

pc pcu^ 

where is the acoustic impedance. Equation (10) can also be expressed 
in dimensionless form as 

= - i 277-77 — P- at y = 0 and 1 (11) 

3p 

3y 


(10) 


It is convenient to express the reciprocal of the impedance ratio in 
terms of the acoustic conductance ratio Kj. and the acoustic susceptance 
ratio a t , that is, 


1 



icr. 


( 12 ) 

( 
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where 1 /£. is called the acoustic admittance ratio 

Finally, substituting Eq. (12) into Eq. (11), and expressing p in 
terms of its real and imaginary terms yields 


ap 


( 1 ) 


ay 


= -27777 


,(U. 


>(2) 


at y = 0 and 1 


(13) 


ao (2) 

— = — 2 7777 

ay 


.( 2 ) 


a + p 


( 1 ) 


K. 


at y = 0 and 1 


(14) 


These are the boundary conditions on pressure to be used at the upper and 
lower surfaces of the duct. 

The impedance at the exit plane, x = L/H, leads to a similar expres- 
sion except an e subscript would be used in place of t in Eqs. (13) and 
(14) and the pressure gradient would be with respect to the variable x. 

Entrance Conditions 

For the entrance pressure profile, the assumption used by Rice 
(ref. 1) of a uniform profile will be used herein. At the present time, 
there is insufficient information available to improve on this assumption. 
For a uniform pressure profile, as will be used in the example considered 
later, 

(0, y) = 1 p^ (0, y) = 0 


(15) 



Axial Acoustic Power 


The sound power which leaves a duct and reaches the far field is 
related to the axial intensity. The axial intensity can be expressed in 
dimensionless form as 

I = i* - —1_ Re (p* u ) (16) 

2 x 

P a 2; r?7 

For hard wall duct with a pc exit impedance, I is identical to 1 for 
all frequencies and duct lengths. 

In terms of the complex representation, the expression for I becomes 


I = 


,(2) . „(1) 3P (2) 1 


2 ITT] 


3x 


9x 


( 17 ) 


The total dimensionless acoustic power is the integral of the intensity 
across the test section 


E x = P I(x, y) dy 


(18) 


The minimum value of E will occur at the duct exit, x = L/H. At this 

A 

position, the sound attenuation will be a maximum. 

By definition, the sound attenuation (the decrease in decibels of the 
acoustic power from x = 0 to x) can be written as 
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Finite Difference Formulation 

Instead of a continuous solution for pressure, the pressure will be 
determined at isolated grid points, as shown in Fig. 2 by means of the 
finite difference approximations. The governing equations and boundary 
conditions can be approximated in difference form (ref. 8) by using either 
a Taylor series expansion, a variational, or an integral formulation. In 
this acoustic problem, where the gradient is specified along a boundary, 
the integration method for generating the finite difference approximations 
is most convenient. 

Wave equation . - The wave equation in finite difference form is 
developed by applying the integration method (ref. 8, p. 168) to the cell 
marked No. 1 in Fig. 2. The details necessary for the development of 
the difference equations are given in Ref. 8 as well as many other texts; 
consequently, only the final results are now presented. The difference 
form of the p^ wave equation, Eq. (6), becomes 



which is the basic Helmholtz equation in difference form. Using the same 

( 2 ) 

procedure, the difference equation for the p' ’ component can also be 
found. 
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Similarity, the difference equation derived from Eq. (13) which 
applies in cell No. 2 in Fig. 1 can be expressed as: 



Similar difference equations exist at the exit cell mark No. 3 and No. 4. 

Axial intensity . - In terms of the difference notation, the axial 
intensity as given by Eq. (17) can be expressed as 


*ii 


<2>L(1) 


2nrj Ax 


i, 3 V 1 , 3 


v i-i, j 


- p! 1 ? 
P i,3 


,(2) _ J2) 


1,3 


• ~ P- i 

1 K i-1 


,3 


(22) 


The total intensity across the test section, as given by Eq. (18), is written 
in difference notation as 


m-1 


E. = 

l 



I. . + - I. 


i,3 


2 l > m 


Ay 


(23) 


By evaluating E. at the entrance and exit positions, taking the log of their 
ratio and multiplying by 10, as indicated in Eq. (19), the maximum sound 
attenuation for the duct is determined. 

We shall now apply the difference equations to the problem of noise 


attenuation in a duct. 
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MATRIX SOLUTION 

The collection of the various difference equations at each grid point 
forms a set of simultaneous equations which can be expressed in matrix 
notation as 

A * P = F (24) 

where A is the coefficient matrix, P is the pressure matrix containing 
the unknown pressures, and F the known matrix containing the various 
initial conditions, where A, P, and F are complex in general. 

The solution to the matrix Eq. (24) yields the values of the real and 
imaginary pressures P^ and P^ at each point. Once these pressures 
are known, the intensity at any position can be determined from Eq. (22), 
and the total intensity can be evaluated from Eq. (23). Finally, the sound 
attenuation in the duct can be found by substituting the values of the total 
intensity into Eq. (19). 

In considering solutions to Eq. (24), it is convenient to express Eq. (24) 
in terms of all real quantities. To accomplish this, the matrix P is 
written as a column vector in terms of the p^ and p^ pressures and 
subdivided as follows: 


A t -C 

1 


'p(D" 


_ F»r 

1 

C | A 


p(2) 

= 

f (2) 

1 J 


_ -*■ 


— — 


( 25 ) 
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The A^ matrix has a form typical of those matrices found in two- 
dimensional heat conduction problems, while the C matrix represents 
the coupling that occurs between the p^ and p^ pressures through 
the impedance boundary conditions. The C matrix is a sparse matrix 

with only one main diagonal term. 

2 

The frequency term (27r?]Ax) in Eq. (20) subtracts from the term 
which will occupy the main diagonal element of the coefficient matrix; 
consequently, for sufficiently high frequencies or spacing parameters, the 
matrix will no longer be positive definite (ref. 8, p. 68). As a result, 
conventional iteration techniques cannot be used. However, matrices of 
the form of Eq. (20) can be solved by elimination techniques. In particular, 
the Gauss elimination technique will be used to find a solution in the 
example problem which now follows. 

Because iteration techniques cannot be used at present, the finite 
difference technique cannot currently be applied to complicated problems 
which require a large number of grid points. Iteration technique or more 
efficient closed form solution will have to be developed to overcome the 
present grid size limitations. 

EXAMPLE CALCULATION 

0 Two example calculations are present in this section so that a com- 
parison between the finite difference solutions and the analytical solutions 
can be made. First, the simple case of plane waves propagating down a 
hard wall duct is presented. This case will allow the comparison of the 
numerical and analytical pressure profiles down the length of the duct. In 
addition, a rule of thumb on the required grid space is developed. The 
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second example will compare the numerical and analytical predictions of 
the attenuation in a soft walled two-dimensional duct. 

ONE- DIMENSIONAL HARD WALL DUCT 

Numerical and analytical values of the pressure were computed for 
the special case of a hard wall duct with a pc exit impedance. Because 
there is no variation of pressure in the y direction, the two-dimensional 
grid lattice shown in Fig. 2 can be reduced to a one -dimensional lattice as 
shown in the upper sketch of Fig. 3. The calculation was made for a hard 
wall (Z = °°) duct with an L/H of 1 and an inlet plane wave with a dimen- 
sionless frequency i] equal to 1. The analytical and numerical values of 
the acoustic pressure profiles along the duct are shown in Fig. 3. As seen 
in Fig. 3, the agreement between the numerical and analytical results is 
good. 

By a series of numerical calculations, the number of grid points 
necessary to get pressure profiles, velocities, and intensity accurate to 
about four percent was found to be 

~ T 

n > 12 77 — 

H 

If n is less than one half the value given by this equation the accuracy 
decreases to 10 percent. The accuracy could be further increased by 
using a four part derivative approximation for the pressure gradient in 
Eq. (17). 
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TWO-DIMENSIONAL SOFT WALLED DUCT 

As another example of the finite difference formulation, the maximum 
noise attenuation will now be calculated for a two-dimensional duct as 
shown in Fig. 2 with L/H of 0. 5 and input plane waves with dimensionless 
frequencies 77 = Hf/c equal to 0. 6 and 1. At present, for two-dimensional 
solutions, the technique is limited to L/H values of 0. 5 or less because 
of matrix size limitations discussed previously. 

The numerical results can be compared to the corresponding analytical 
results found by using the techniques presented by Rice (ref. 1), which is 
applicable to infinite ducts (no reflected waves). In the numerical formula- 
tion of the infinite duct, the duct can be simulated in a finite length by 
using an impedance value of pc (k 0 = 1; = 0) at the exit plane of the 

duct. For plane wave propagating in a hard-wall duct, a pc exit impedance 
leads to zero reflected wave which is analogous to an infinitely long duct. 
This assumed exit condition is discussed in greater detail in the Appendix 
of this paper. 

The first step in determining the maximum attenuation in the duct is 
to calculate the attenuation at various Z values in the impedance plane. 

This is shown in Fig. 4 where 

C = — = 0+ i X (26) 

pc 

The sound attenuation values shown in Fig. 4 were obtained by choosing 
discrete values of 8 and x throughout the 6-\ plane, calculating the 
appropriate a and k from Eq. (12), solving the matrix Eq. (25) by 
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Gauss elimination for the pressure distribution, and finally calculating the 
sound attenuation at the duct exit from Eqs. (23) and (19). The results 
shown in Fig. 4 were obtained by interpolating between the points to obtain 
smooth contours. The impedance associated with the peak attenuation is 
(.called the optimum impedance. 

1 

Because of matrix size limitations in the subprogram used to solve 
Eq. (25), the technique of varying grid size (ref. 9) was used. This is 
illustrated in Fig. 5 where the grid spacing Ay is decreased toward zero. 
The attenuation at the optimum point is found by extrapolating the grid size 
to zero (dotted portion of the curve). As seen in Fig. 5, the extrapolated 
numerical value is in agreement with the analytical values calculated from 
the theory presented by Rice in Ref. 1. • 

This agreement is very encouraging considering the few grid points 
(maximum of 50) used in obtaining the results. Once new matrix iteration 
or elimination techniques are developed, the prognosis for obtaining accurate 
results for complex problems with a minimum of computer storage and 
running time is excellent. 

STEADY FLOW PROBLEM 

At present, I am extending the difference technique to the steady flow 
case. Because of the added Mach number terms in both the wave equation 
and the slip boundary condition, the coupling matric C in Eq. (25) becomes 
tridiagonal. However, these complications have not as yet posed any 
barriers in obtaining the steady flow solutions. 


4 ' 
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CONCLUDING REMARKS 

A finite difference solution for sound propagation in a two-dimensional 
soft-walled duct for zero Mach number shows favorable comparison with 
the corresponding exact analytical results. Because the solution matrix 
for the acoustical flow field is not positive definite, conventional iteration 
techniques cannot be used to solve the difference equations. Before the 
finite differences formulation can be applied to more complicated problems 
which require a large number of grid points, such as the sonic inlet, 
iteration techniques will have to be developed to overcome present grid 
size limitations. 

The finite difference formulation is flexible and should be a powerful 
tool in the solution of more realistic studies of inlet and exhaust ducts of 
turbofan engines. The present formulation allows complete freedom in 
choosing the inlet pressure profile and the complex impedance along both 
boundaries. The extension of the present formulation to both uniform flow 
and shear flow is straight forward. 



APPENDIX - EXIT IMPEDANCE ASSUMPTION 


The close agreement between the numerical theory and the correspond- 
ing analytical results substantiates the utility of the pc exit impedance 
assumption made in the previous section of this report. However, the 
analytical theories (ref. 1 for example) indicate that certain alterations in 
the numerical exit conditions are required if the numerical theory is to be 
in exact agreement with the analysis. 

The analytical theory for an infinite duct with no reflected waves 
assumes the solutions to the wave equation are separable and can be 
expressed in terms of an infinite number of modes. However, the infinite 
duct with no reflections could be replaced by a finite length of duct L with 
the exit impedance at L chosen so that no reflections would occur. For 
no reflections, the analysis indicates that each mode has its own unique 
exit impedance. Thus, reflections must be occurring at the duct exit in 
the numerical calculation since choosing an exit impedance equal to pc 
can not possibly simulate all the distinct impedance values required for 
each mode in the analysis. 

Fortunately, in many practical problems, most of the higher order 
modes decay out and the lowest order mode appears at the exit. A com- 
parison to the analytical results indicates that, in general, the exit 
impedance associated with the lowest order mode is close to pc. Never- 
theless, it is important to establish a general numerical procedure which 
will converge to the correct answer. 

Overlooking round-off and truncation errors in the numerical calcu- 
lations, the attenuation calculated by the numerical technique will be greater 
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than that calculated by the exact analytical technique because of reflec- 
tions which occur at the exit; that is 

AdB\ > AdB\ 

Inumerical (analytical 

7o<x< L / 0 < x < L 

The reflected energy appears to have been absorbed by the soft walls in 
the numerical calculation. 

Now, consider the case where an additional length of absorbing 
liner AL is added to the original liner, so that the liner now has a 
length L + AL. If AL is chosen sufficiently long, the reflected energy 
at the duct exit x = L + AL will be absorbed before this energy can 
reach the previous exit position at x = L. Therefore, for soft-walled 
ducts 


limit AdB 
AL — 00 


numerical 


0 < x < L 


AdB\ 


| analytical 
<x< L 


I (U 

f 0 


Practically, the liner length need only be increased by some minimum 
length AL until the attenuation for the liner length between x = 0 and 
x = L remains constant to a given percentage. 
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Figure Z - Coordinate and grid point representation of two dimensional soft walled duct 
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Figure 3. - Analytical and numerical pressure profiles for one 
dimensional sound propagation with pc exit impedance for 
i? » 1 with n ■ 12. 



Figure 4 - Sound power attenuation contours for rj • 
0. 5 and m - 20, n - 5. 
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Figure 5. - Effect of spacing on attenuation at optimum im- . 
pedance in two dimensional duct for L/H - 0. 5 where 
7) 3 0. 6 with n 3 5 and r\ = 1 with n 3 10. 
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